#########################################################################
#  File-Name:	cyprus-deactivation.R					
#  Date: September 2023
#  Purpose: This file analyzes (survey) results from the deactivation RCT in Cyprus
# 
#  Data In: cyprus-replication-data-2023.csv
#  Data Out: cyprus-replication-data-2023.csv
#  Status: Completed						
#  Machine:  Nejla Asimovic Local						
#########################################################################

library(ggplot2); library(ggpubr); 
library(stringr); library(estimatr); 
library(Rmisc); library(devtools);
library(plotly); library(dplyr); 
library(sandwich); library(plyr);
library(Rmisc); library(lmtest); 
library(tibble);  library(stringi)
library(estimatr); library(xtable)
library(knitr); library(stringr)
library(kableExtra); library(expss)
library(rio); library(car);
library(ggplot2); library(arm)

# Set your own working directory
rootArg <-  if(Sys.info()["user"]=="nejlaasimovic"){
  setwd("~/Desktop/desktop-sept23/final-r-p")
}

# Upload dataset from wherever you have it saved
data <-  read.csv("cyprus-replication-data-2023.csv") 
data <- data[data$endline_survey=="1",] # limiting analysis to those whose valid final survey responses we obtained
options(digits=3)

# Source the analysis functions
source("~/Desktop/desktop-sept23/final-r-p/functions-cyprus-2023.R")

# -----------------------------------------
#   MAIN RESULTS: INTER-ETHNIC ATTITUDES
# -----------------------------------------

# Imputing modes for missing variables at a gender-ethnicity level, as in pre-registration

data$gender_ethn <- paste0(data$gender,
                           "-",
                           data$ethn)
covariates <- c("fb_weeklyusage","gender","age","educ","imp_ethn","imp_cntry")
imputing <- cbind(sapply(covariates, function (x)  imputeGroupMode(data, x,"gender_ethn")))


# ------- AFFECT/FEELING THERMOMETER

# Creating feeling thermometer variable: subtracting ft of the outgroup from ft of the ingroup
data$ft <- data$outgroup_ft-data$ingroup_ft # recoded, so that higher number refers to a more positive atribution

# ------- COOPERATION 
# Creating cooperation index as a sum of z-scales from the disagreement with three statements:
# 1.Multi-ethnic parties cannot secure the interest of my people.
# 2.It is only possible to cooperate with members of my ethnic group.

data$cooperation <- NA
for (i in 1:nrow(data)){
  data$cooperation[i] <-data$cooperation_1[i]+ data$cooperation_2[i]
}


#  -------- PERCEPTION OF OUTGROUP EVALUATIONS
# higher values indicate higher level of agreement, hence negative characteristics need to be recoded (multiplied by -1 within the index creation)

data$inchrct_data<-NA
for (i in 1:nrow(data)){
  data$inchrct_data[i] <- data$in_broad[i] - data$in_hyp[i]- data$in_evil[i] + data$in_patriot[i] + data$in_gener[i]  + data$in_honest[i] - data$in_self[i] - data$in_rude[i]}

#  -------- OUTGROUP CHARACTERISTICS
# higher values indicate higher level of agreement, hence negative characteristics need to be recoded (multiplied by -1 within the index creation)

data$othchrct_data<-NA
for (i in 1:nrow(data)){
  data$othchrct_data[i] <- data$oth_broadv[i] - data$oth_hyp[i] - data$oth_evil[i] + data$oth_patriot[i] + data$oth_gener[i]  + data$oth_honest[i] - data$oth_self[i] - data$oth_rude[i]}

# --------- SOCIAL DISTANCE

data$sd <- NA
for (i in  1:nrow(data)){
  data$sd[i] <- data$sd_1_marry[i]  + data$sd_2_closefriends[i] +data$sd_3_neighbor[i] + data$sd_4_colleague[i] + data$sd_5_ascitizens[i]}

# ------ OUTGROUP INDEX
# -- Creating the index as a sum of z-scores
data$outgroup_index  <- NA
data$outgroup_index <-  scale(data$sd) + scale(data$cooperation) + scale(data$ft) + scale(data$othchrct_data)+scale(data$inchrct_data)

# -- Creating the index as a PC index

outgroup_variables <- c("sd",
                        "cooperation",
                        "inchrct_data",
                        "ft","othchrct_data")
outgroup_variables_pc <- data[,outgroup_variables] 
psych::alpha(outgroup_variables_pc, 'check.keys=TRUE')$total$std.alpha # checking internal validity and Chronbach alpha
data$outgroup_princomp <- pc_index(dataIn=data, varList=outgroup_variables)


# ------------- FIGURE 1 -----------------

coef.vec1<- c(reg_res_mod3("ft",data)["estimate3"], reg_res_mod3("sd",data)["estimate3"],reg_res_mod3("cooperation",data)["estimate3"],
              reg_res_mod3("inchrct_data", data)["estimate3"], reg_res_mod3("othchrct_data", data)["estimate3"], reg_res_mod3("outgroup_index",data)["estimate3"])
se.vec1<- c(reg_res_mod3("ft",data)["std.error3"], reg_res_mod3("sd",data)["std.error3"],reg_res_mod3("cooperation",data)["std.error3"],
            reg_res_mod3("inchrct_data", data)["std.error3"], reg_res_mod3("othchrct_data", data)["std.error3"], reg_res_mod3("outgroup_index",data)["std.error3"])

coef.vec2<- c(reg_res_mod3y("ft",data[data$year=="20",])["estimate3"], reg_res_mod3y("sd",data[data$year=="20",])["estimate3"],reg_res_mod3y("cooperation",data[data$year=="20",])["estimate3"],
              reg_res_mod3y("inchrct_data",data[data$year=="20",])["estimate3"], reg_res_mod3y("othchrct_data",data[data$year=="20",])["estimate3"], reg_res_mod3y("outgroup_index",data[data$year=="20",])["estimate3"])
se.vec2<- c(reg_res_mod3y("ft",data[data$year=="20",])["std.error3"], reg_res_mod3y("sd",data[data$year=="20",])["std.error3"],reg_res_mod3y("cooperation",data[data$year=="20",])["std.error3"],
            reg_res_mod3y("inchrct_data",data[data$year=="20",])["std.error3"], reg_res_mod3y("othchrct_data",data[data$year=="20",])["std.error3"], reg_res_mod3y("outgroup_index",data[data$year=="20",])["std.error3"])

coef.vec3<- c(reg_res_mod3y("ft",data[data$year=="21",])["estimate3"], reg_res_mod3y("sd",data[data$year=="21",])["estimate3"],reg_res_mod3y("cooperation",data[data$year=="21",])["estimate3"],
              reg_res_mod3y("inchrct_data",data[data$year=="21",])["estimate3"], reg_res_mod3y("othchrct_data",data[data$year=="21",])["estimate3"], reg_res_mod3y("outgroup_index",data[data$year=="21",])["estimate3"])
se.vec3<- c(reg_res_mod3y("ft",data[data$year=="21",])["std.error3"], reg_res_mod3y("sd",data[data$year=="21",])["std.error3"],reg_res_mod3y("cooperation",data[data$year=="21",])["std.error3"],
            reg_res_mod3y("inchrct_data",data[data$year=="21",])["std.error3"], reg_res_mod3y("othchrct_data",data[data$year=="21",])["std.error3"], reg_res_mod3y("outgroup_index",data[data$year=="21",])["std.error3"])

dev.off()
adjust = 0.10
layout(matrix(c(2,1),1,2.5), #in order to add variable categories and braces to left side of plot, 
       widths = c(2, 5.8))
var.names <- c("Feeling thermometer","Social distance","Cooperation","Metastereotypes","Outgroup ratings","OUTGROUP INDEX")
y.axis <- c(length(var.names):1) 
plot(coef.vec1, y.axis, type = "p", axes = F, xlab = "Standardized effect of deactivation", ylab = "", pch = 19,cex = 1,#plot coefficients as points, turning off axes and labels.
     xlim = c(-1,1), xaxs = "r", main = "Outgroup regard", col="red", cex.lab=1.1,mgp = c(3,1,0.2))
segments(unlist(as.numeric(coef.vec1))-qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, unlist(as.numeric(coef.vec1))+qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, lwd = 1, col = c("red"))
axis(1, at = seq(-1,1,by=0.5), labels =c(-1,-0.5, 0,0.5,1), tick = T,#draw x-axis and labels with tick marks
     cex.axis = 1, mgp = c(4,2,1))#reduce label size, moves labels closer to tick marks; zadnji u mgp spusta x liniju nize
axis(2, at = y.axis, label = var.names, las = 1, tick = T, mgp = c(1.6,.7,0),cex.axis = 1.35) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
abline(h = 0, v=0, col="gray",lty=2,lwd=2.5)
segments(unlist(as.numeric(coef.vec2))-qnorm(.975)*unlist(as.numeric(se.vec2)), y.axis-adjust, unlist(as.numeric(coef.vec2))+qnorm(.975)*unlist(as.numeric(se.vec2)), y.axis-adjust,lwd = 1, col = c("darkgrey"))#draw lines connecting 95% confidence intervals
points(coef.vec2, y.axis-adjust,pch = 17, cex = 0.8, col = c("darkgrey")) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

adjust=0.15
segments(unlist(as.numeric(coef.vec3))-qnorm(.975)*unlist(as.numeric(se.vec3)), y.axis-adjust, unlist(as.numeric(coef.vec3))+qnorm(.975)*unlist(as.numeric(se.vec3)), y.axis-adjust,lwd = 1, col = c("black"))#draw lines connecting 95% confidence intervals
points(coef.vec3, y.axis-adjust,pch = 23, cex = 0.8, col = c("black")) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color
legend(x=0.5, y=6.5,
       legend = c("Pooled", "2020 Wave","2021 Wave"),
       col = c("red", "darkgrey","black"), xpd=TRUE,
       pch = c(19,17, 23),
       bty = T,
       pt.cex = 1,
       cex = 0.8,
       text.col = "black",   box.lty=1, box.lwd=1,
       horiz=F)

# -----------------------------------------------
#       Facebook Substitutes -  Fig 2 
# -----------------------------------------------
# Left figure
dev.off()
par(mar=c(8,12,5,5))
coef.vec<- c(reg_res_mod2("activ_1", data)["estimate2"],
             reg_res_mod2("activ_2", data)["estimate2"],
             reg_res_mod2("activ_3", data)["estimate2"],
             reg_res_mod2("activ_4", data)["estimate2"],
             reg_res_mod2("activ_5", data)["estimate2"])
se.vec<- c(reg_res_mod2("activ_1", data)["std.error2"],
           reg_res_mod2("activ_2", data)["std.error2"],
           reg_res_mod2("activ_3", data)["std.error2"],
           reg_res_mod2("activ_4", data)["std.error2"],
           reg_res_mod2("activ_5", data)["std.error2"])
var.names <- c("Instagram","Other online activities","TV","Family & friends","No technology")
results_df <- data.frame(term=var.names, estimate=coef.vec,
                         std.error=se.vec)
y.axis <- c(length(coef.vec):1) 
plot(coef.vec, y.axis, type = "p", axes = F, xlab = "Treatment effect (SD)", cex.lab=1.35,ylab = "", pch = 19,cex = 1.5,
     xlim = c(-1,1), xaxs = "r", main = "Facebook substitutes", cex.main=1.35, mgp = c(3,.7,0.2), col=c("black",
                                                                                                        "black",
                                                                                                        "black",
                                                                                                        "black",
                                                                                                        "black" ),cex.lab=1.15)
axis(1, at = seq(-1,1,by=0.5), labels =c(-1,-0.5,0,0.5,1) , tick = T,#draw x-axis and labels with tick marks
     cex.axis = 1.2, mgp = c(4,2,1))
axis(2, at = y.axis, label = var.names, las = 1.2, tick = T, mgp = c(1.6,.7,0),cex.axis = 1.4) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
abline(h = 0, v=0, col="gray",lty=2,lwd=2.5)
segments(unlist(coef.vec)-qnorm(.975)*unlist(se.vec), y.axis, unlist(coef.vec)+qnorm(.975)*unlist(se.vec), y.axis,lwd = 1.5, col = c("black"))

# ---- Right figure (barplot of activities when deactivated)
inst <- mean(data[data$treat=="1",]$activ_1, na.rm=T)
online <- mean(data[data$treat=="1",]$activ_2, na.rm=T)
tv <- mean(data[data$treat=="1",]$activ_3, na.rm=T)
fam <- mean(data[data$treat=="1",]$activ_4, na.rm=T)
offscreen <- mean(data[data$treat=="1",]$activ_5, na.rm=T)
pal <- colorRampPalette(colors = c("purple", "blue"))(7)

dev.off()
# png("figure2.png",   # The directory you want to save the file in
# width = 7, # The width of the plot in inches
# height = 4.1) 
par(mar = c(10, 30, 5, 10))
x <- c("Friends & Family", "Instagram","No Tech", "Other online","TV")
y <- c(fam, inst,offscreen,online,tv)
data_x <- data.frame(x, y, stringsAsFactors = FALSE)
data_x$x <- factor(data_x$x, levels = unique(data_x$x)[order(data_x$y, decreasing = TRUE)])
plot_ly(data_x,x=~x, y=~y, 
        name="FDs",type="bar",
        marker = list(color = c("rgba(70, 120, 150, 5)",
                                "rgba(70, 120, 150, 0.4)",
                                "rgba(70, 120, 150, 0.9)",
                                "rgba(50, 120, 150, 0.7)",
                                "rgba(50, 120, 150, 0.2)",
                                "rgba(50, 120, 150, 0.2"))) %>% layout(xaxis = list(title = "", tickangle = -45,categoryorder = "array",tickfont=list(size=19)),
                                                                       yaxis = list(title = "Average response",categoryorder = "array",tickfont=list(size=19), mgp=c(2,1,0),titlefont = list(size=20),line=2, ylim(1,5)))



# ------------------------
#    OFFLINE NETWORKS 
# ------------------------
# offline_network - outgroup friends
# offline_exposure - passive exposure
# offline_active - active exposure

data$offline_network <- 6-data$offline_network # reorient so that higher value indicates more homogenous network
data$offline_exposure <- 6-data$offline_exposure
data$offline_active <- 6-data$offline_active

# Creating offline network index as a sum score
data$offline_network_diversity <- NA
data[data$year=="21",]$offline_network_diversity <- scale(data[data$year=="21",]$offline_network) + scale(data[data$year=="21",]$offline_exposure) + scale(data[data$year=="21",]$offline_active)

# Creating offline network index as a PC score
network_variables <- c("offline_network","offline_exposure","offline_active")
network_variables_pc <- data[,network_variables]
psych::alpha(network_variables_pc, 'check.keys=TRUE')$total$std.alpha 
network_variables_pc <- data[data$year=="21",][,network_variables]
data$offline_network_diversity_pc<- pc_index(dataIn=data, varList=network_variables)

# --------------------- Table 2 --------------
# Column 1: Offline homogenous network (sum score)
interaction <- lm(outgroup_index ~ treat + offline_network_diversity + treat:offline_network_diversity, data[data$year=="21",])

interaction_est_coef <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[2,1]
interaction_est_se <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[2,2]
interaction_est_p <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[2,4]

interaction_int_coef <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[4,1]
interaction_int_se <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[4,2]
interaction_int_p <- coeftest(interaction, vcov = vcovHC(interaction, type = "HC1"))[4,4]

# Column 2: Offline homogenous network (PC)
interaction_pc <- lm(outgroup_index ~ treat + offline_network_diversity_pc   + treat:offline_network_diversity_pc, data[data$year=="21",])

interaction_pc_est_coef <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[2,1]
interaction_pc_est_se <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[2,2]
interaction_pc_est_p <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[2,4]

interaction_pc_int_coef <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[4,1]
interaction_pc_int_se <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[4,2]
interaction_pc_int_p <- coeftest(interaction_pc, vcov = vcovHC(interaction_pc, type = "HC1"))[4,4]

# Putting it together
coef.vec_treat <- c(interaction_est_coef, interaction_pc_est_coef)
coef.vec_intr <- c(interaction_int_coef, interaction_pc_int_coef)

se.vec_treat <- c(interaction_est_se,interaction_pc_est_se)
se.vec_intr <- c(interaction_int_se, interaction_pc_int_se)

p.vec_treat <-  c(interaction_est_p, interaction_pc_est_p)
p.vec_intr <- c(interaction_int_p,interaction_pc_int_p)

names <- c("Outgroup regard (sum score)","Outgroup regard (PC)")
merged_intr <- data.frame(names, coef.vec_treat,se.vec_treat, p.vec_treat,coef.vec_intr, se.vec_intr, p.vec_intr)
merged_intr  # Table 2

# ----------------------------
#           APPENDIX
# ----------------------------

# --------------- SAMPLE DESCRIPTION ---------------
# ---------------- Table A1 ------------------------

table(data$year)
table(data$treat)
table(data[data$year=="20",]$treat)
table(data[data$year=="21",]$treat)

# ------------- Table A2 ------------------
dat_cyp <- data[!is.na(data$cypriots) & data$lives_not_cyprus=="0",] # 442
dat_diasp <- data[!is.na(data$cypriots) & data$lives_not_cyprus=="1",] # 61
dat_immigrants <- data[is.na(data$cypriots) & is.na(data$lives_not_cyprus),] # [taken out as explained in the Appedix; previously 22]
dat_homelands <- data[is.na(data$cypriots) & !is.na(data$lives_not_cyprus),] # 40

# ----------------------------------------
#       TREATMENT vs. CONTROL
# ----------------------------------------
# ------------- Table A3 ------------------

# Here we compare baseline characteristics between treatment and control group:
bal.variables <- c("age","gender","educ","employ","fbuser_length","fb_weeklyusage","fb_dailyusage","news_freq", "pol_int", "outgroup_contact", "numb_forums","trust_media", "imp_ethn","imp_cntry","imp_relg")
df <- data[,bal.variables]

p.values <- cbind(sapply(bal.variables, function (x) t.test(data[data$treat==1,][,x], data[data$treat==0,][,x]))[3,])
balance_table <- do.call(data.frame, # creating a table combining the mean values and p-values
                         list(TreatmentGroup = apply(df[data$treat == "1",], 2, mean, na.rm=TRUE),
                              ControlGroup =  apply(df[data$treat =="0",], 2, mean, na.rm=TRUE),
                              TreatmentSD =  apply(df[data$treat == "1",], 2, sd, na.rm=TRUE),
                              ControlSD =  apply(df[data$treat =="0" ,], 2, sd, na.rm=TRUE),
                              t.test=p.values))
balance_table 

# ----------------------------------------
#              ATTRITION 
# ----------------------------------------
attrition <- read.csv("cyprus-replication-attrition-2023.csv")
table(attrition$treat)

for(varUp in bal.variables){
  attrition[,paste0(varUp)] <- as.numeric(attrition[,paste0(varUp)])
}

# ------------- Table A4 ------------------
# attriters: treatment vs. control
p.values <- cbind(sapply(bal.variables, function (x) t.test(df[attrition$treat==1,][,x], df[attrition$treat==0,][,x]))[3,])
balance_table <- do.call(data.frame, # creating a table combining the mean values and p-values
                         list(TreatmentGroup = apply(df[attrition$treat == "1",], 2, mean, na.rm=TRUE),
                              ControlGroup =  apply(df[attrition$treat =="0",], 2, mean, na.rm=TRUE),
                              TreatmentSD =  apply(df[attrition$treat == "1",], 2, sd, na.rm=TRUE),
                              ControlSD =  apply(df[attrition$treat =="0" ,], 2, sd, na.rm=TRUE),
                              t.test=p.values))
balance_table

# ------------- Table A5 ------------------
# final data vs. attriters
p.values <- cbind(sapply(bal.variables, function (x) t.test(data[,x], attrition[,x]))[3,])
balance_table <- do.call(data.frame, # creating a table combining the mean values and p-values
                         list(Final = apply(data[,bal.variables], 2, mean, na.rm=TRUE),
                              Attrit =  apply(attrition[,bal.variables], 2, mean, na.rm=TRUE),
                              FinalSD =  apply(data[,bal.variables], 2, sd, na.rm=TRUE),
                              AttritSD =  apply(attrition[,bal.variables], 2, sd, na.rm=TRUE),
                              t.test=p.values))
balance_table

# Does treatment status predict whether a user finished the survey or not?
summary(lm(endline_survey~treat, data=data)) # Section 2 | "Attrition"


# ------------- RESULTS SECTION ------------------
# For replicating the tables within the results section, we first need to calculate
# treatment effects on subjective well-being and news knowledge.

# -------- Subjective WellBeing
# Variables are coded so that higher values indicate higher level of agreement;
# hence, negative variables need to be recoded

data$nerv <- (-1) * data$nerv 
data$depression <- (-1) * data$depression
data$loneliness <- (-1) * data$loneliness
data$boredom <- (-1) * data$boredom
data$isol<- (-1) * data$isol

for (i in 1:nrow(data)){
  data$swb[i] <- scale(data$satisf)[i]+scale(data$joy)[i]+scale(data$fulf)[i]+scale(data$depression)[i] +scale(data$loneliness)[i] + scale(data$nerv)[i] 
  +scale(data$boredom)[i]+scale(data$isol[i])
}



# -------- News knowledge
c_news<-NA
for (i in 1:nrow(data)){
  data$c_news[i] <- data$fbnews_st1[i]+data$fbnews_st2[i]+data$fbnews_st3[i]+data$fbnews_st4[i]+data$fbnews_st5[i]+data$fbnews_st6[i]+data$fbnews_st7[i]+data$fbnews_st8[i]
}

# ---------------------------
#     TABLE A7
# ---------------------------

# ----- Top table A7: news knowledge + subjective well-being
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1("c_news",data),
  reg_res_mod1("swb",data),
  reg_res_mod1("satisf",data),
  reg_res_mod1("joy",data),
  reg_res_mod1("fulf",data),
  reg_res_mod1("nerv",data),
  reg_res_mod1("boredom",data),
  reg_res_mod1("loneliness",data),
  reg_res_mod1("depression",data),
  reg_res_mod1("isol",data)))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2("c_news",data),
  reg_res_mod2("swb",data),
  reg_res_mod2("satisf",data),
  reg_res_mod2("joy",data),
  reg_res_mod2("fulf",data),
  reg_res_mod2("nerv",data),
  reg_res_mod2("boredom",data),
  reg_res_mod2("loneliness",data),
  reg_res_mod2("depression",data),
  reg_res_mod2("isol",data)))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3("c_news",data),
  reg_res_mod3("swb",data),
  reg_res_mod3("satisf",data),
  reg_res_mod3("joy",data),
  reg_res_mod3("fulf",data),
  reg_res_mod3("nerv",data),
  reg_res_mod3("boredom",data),
  reg_res_mod3("loneliness",data),
  reg_res_mod3("depression",data),
  reg_res_mod3("isol",data)))
Table_MainResults3

tableA7_top <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
tableA7_top # Top Table A7: news knowledge and subjective well-being

# ----- Bottom table A7 : outgroup regard
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1("sd",data),
  reg_res_mod1("cooperation",data),
  reg_res_mod1("ft",data),
  reg_res_mod1("othchrct_data",data),
  reg_res_mod1("inchrct_data",data),
  reg_res_mod1("outgroup_index",data),
  reg_res_mod1("outgroup_princomp",data)))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2("sd",data),
  reg_res_mod2("cooperation",data),
  reg_res_mod2("ft",data),
  reg_res_mod2("othchrct_data",data),
  reg_res_mod2("inchrct_data",data),
  reg_res_mod2("outgroup_index",data),
  reg_res_mod2("outgroup_princomp",data)))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3("sd",data),
  reg_res_mod3("cooperation",data),
  reg_res_mod3("ft",data),
  reg_res_mod3("othchrct_data",data),
  reg_res_mod3("inchrct_data",data),
  reg_res_mod3("outgroup_index",data),
  reg_res_mod3("outgroup_princomp",data)))
Table_MainResults3

TableA7_bottom <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
TableA7_bottom # Bottom Table A7: outgroup regard

# ---------------------------
#     TABLE A8
# ---------------------------

# ----- Top table A8 [2020]: news knowledge + subjective well-being
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1y("c_news",data[data$year=="20",]),
  reg_res_mod1y("swb",data[data$year=="20",]),
  reg_res_mod1y("satisf",data[data$year=="20",]),
  reg_res_mod1y("joy",data[data$year=="20",]),
  reg_res_mod1y("fulf",data[data$year=="20",]),
  reg_res_mod1y("nerv",data[data$year=="20",]),
  reg_res_mod1y("boredom",data[data$year=="20",]),
  reg_res_mod1y("loneliness",data[data$year=="20",]),
  reg_res_mod1y("depression",data[data$year=="20",]),
  reg_res_mod1y("isol",data[data$year=="20",])))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2y("c_news",data[data$year=="20",]),
  reg_res_mod2y("swb",data[data$year=="20",]),
  reg_res_mod2y("satisf",data[data$year=="20",]),
  reg_res_mod2y("joy",data[data$year=="20",]),
  reg_res_mod2y("fulf",data[data$year=="20",]),
  reg_res_mod2y("nerv",data[data$year=="20",]),
  reg_res_mod2y("boredom",data[data$year=="20",]),
  reg_res_mod2y("loneliness",data[data$year=="20",]),
  reg_res_mod2y("depression",data[data$year=="20",]),
  reg_res_mod2y("isol",data[data$year=="20",])))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3y("c_news",data[data$year=="20",]),
  reg_res_mod3y("swb",data[data$year=="20",]),
  reg_res_mod3y("satisf",data[data$year=="20",]),
  reg_res_mod3y("joy",data[data$year=="20",]),
  reg_res_mod3y("fulf",data[data$year=="20",]),
  reg_res_mod3y("nerv",data[data$year=="20",]),
  reg_res_mod3y("boredom",data[data$year=="20",]),
  reg_res_mod3y("loneliness",data[data$year=="20",]),
  reg_res_mod3y("depression",data[data$year=="20",]),
  reg_res_mod3y("isol",data[data$year=="20",])))
Table_MainResults3

tableA8_top <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
tableA8_top # Top Table A8 [2020]: news knowledge and subjective well-being

# ----- Bottom table A8 [2020]: outgroup regard
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1y("sd",data[data$year=="20",]),
  reg_res_mod1y("cooperation",data[data$year=="20",]),
  reg_res_mod1y("ft",data[data$year=="20",]),
  reg_res_mod1y("othchrct_data",data[data$year=="20",]),
  reg_res_mod1y("inchrct_data",data[data$year=="20",]),
  reg_res_mod1y("outgroup_index",data[data$year=="20",]),
  reg_res_mod1y("outgroup_princomp",data[data$year=="20",])))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2y("sd",data[data$year=="20",]),
  reg_res_mod2y("cooperation",data[data$year=="20",]),
  reg_res_mod2y("ft",data[data$year=="20",]),
  reg_res_mod2y("othchrct_data",data[data$year=="20",]),
  reg_res_mod2y("inchrct_data",data[data$year=="20",]),
  reg_res_mod2y("outgroup_index",data[data$year=="20",]),
  reg_res_mod2y("outgroup_princomp",data[data$year=="20",])))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3y("sd",data[data$year=="20",]),
  reg_res_mod3y("cooperation",data[data$year=="20",]),
  reg_res_mod3y("ft",data[data$year=="20",]),
  reg_res_mod3y("othchrct_data",data[data$year=="20",]),
  reg_res_mod3y("inchrct_data",data[data$year=="20",]),
  reg_res_mod3y("outgroup_index",data[data$year=="20",]),
  reg_res_mod3y("outgroup_princomp",data[data$year=="20",])))
Table_MainResults3

tableA8_bottom <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
tableA8_bottom # Bottom Table A8 [2020]: outgroup regard

# ---------------------------
#         TABLE A9
# ---------------------------

# ----- Top table A9 [2021]: news knowledge + subjective well-being
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1y("c_news",data[data$year=="21",]),
  reg_res_mod1y("swb",data[data$year=="21",]),
  reg_res_mod1y("satisf",data[data$year=="21",]),
  reg_res_mod1y("joy",data[data$year=="21",]),
  reg_res_mod1y("fulf",data[data$year=="21",]),
  reg_res_mod1y("nerv",data[data$year=="21",]),
  reg_res_mod1y("boredom",data[data$year=="21",]),
  reg_res_mod1y("loneliness",data[data$year=="21",]),
  reg_res_mod1y("depression",data[data$year=="21",]),
  reg_res_mod1y("isol",data[data$year=="21",])))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2y("c_news",data[data$year=="21",]),
  reg_res_mod2y("swb",data[data$year=="21",]),
  reg_res_mod2y("satisf",data[data$year=="21",]),
  reg_res_mod2y("joy",data[data$year=="21",]),
  reg_res_mod2y("fulf",data[data$year=="21",]),
  reg_res_mod2y("nerv",data[data$year=="21",]),
  reg_res_mod2y("boredom",data[data$year=="21",]),
  reg_res_mod2y("loneliness",data[data$year=="21",]),
  reg_res_mod2y("depression",data[data$year=="21",]),
  reg_res_mod2y("isol",data[data$year=="21",])))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3y("c_news",data[data$year=="21",]),
  reg_res_mod3y("swb",data[data$year=="21",]),
  reg_res_mod3y("satisf",data[data$year=="21",]),
  reg_res_mod3y("joy",data[data$year=="21",]),
  reg_res_mod3y("fulf",data[data$year=="21",]),
  reg_res_mod3y("nerv",data[data$year=="21",]),
  reg_res_mod3y("boredom",data[data$year=="21",]),
  reg_res_mod3y("loneliness",data[data$year=="21",]),
  reg_res_mod3y("depression",data[data$year=="21",]),
  reg_res_mod3y("isol",data[data$year=="21",])))
Table_MainResults3

tableA9_top <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
tableA9_top # Top Table A9 [2021]: news knowledge and subjective well-being

# ----- Bottom table A9 [2021]: outgroup regard
Table_MainResults1 <- do.call(data.frame,rbind(
  reg_res_mod1y("sd",data[data$year=="21",]),
  reg_res_mod1y("cooperation",data[data$year=="21",]),
  reg_res_mod1y("ft",data[data$year=="21",]),
  reg_res_mod1y("othchrct_data",data[data$year=="21",]),
  reg_res_mod1y("inchrct_data",data[data$year=="21",]),
  reg_res_mod1y("outgroup_index",data[data$year=="21",]),
  reg_res_mod1y("outgroup_princomp",data[data$year=="21",])))
Table_MainResults1

Table_MainResults2 <- do.call(data.frame,rbind(
  reg_res_mod2y("sd",data[data$year=="21",]),
  reg_res_mod2y("cooperation",data[data$year=="21",]),
  reg_res_mod2y("ft",data[data$year=="21",]),
  reg_res_mod2y("othchrct_data",data[data$year=="21",]),
  reg_res_mod2y("inchrct_data",data[data$year=="21",]),
  reg_res_mod2y("outgroup_index",data[data$year=="21",]),
  reg_res_mod2y("outgroup_princomp",data[data$year=="21",])))
Table_MainResults2

Table_MainResults3 <- do.call(data.frame,rbind(
  reg_res_mod3y("sd",data[data$year=="21",]),
  reg_res_mod3y("cooperation",data[data$year=="21",]),
  reg_res_mod3y("ft",data[data$year=="21",]),
  reg_res_mod3y("othchrct_data",data[data$year=="21",]),
  reg_res_mod3y("inchrct_data",data[data$year=="21",]),
  reg_res_mod3y("outgroup_index",data[data$year=="21",]),
  reg_res_mod3y("outgroup_princomp",data[data$year=="21",])))
Table_MainResults3

TableA9_bottom <-join_all(list(Table_MainResults1, Table_MainResults2, Table_MainResults3), by="outcome")
TableA9_bottom  # Bottom Table A9 [2021]: outgroup regard

# --------------------------------------
#     FDR Adjusted Results
# --------------------------------------

# Multiple comparison adjustment (Benjamini-Hochberg)

# ------------- Table A10 [news knowledge and subjective well-being]------------------
p.values_swb <- as.table(c(
  reg_res_mod3("satisf",data)$p.value,
  reg_res_mod3("joy",data)$p.value,
  reg_res_mod3("fulf",data)$p.value,
  reg_res_mod3("nerv",data)$p.value,
  reg_res_mod3("boredom",data)$p.value,
  reg_res_mod3("loneliness",data)$p.value,
  reg_res_mod3("depression",data)$p.value,
  reg_res_mod3("isol",data)$p.value,
  reg_res_mod3("swb",data)$p.value))

names <- c("Satisfaction","Joy","Fulfillment","Anxiety (rev coded)","Boredom (rev coded)","Loneliness (rev coded)", "Depression (rev coded)","Isolation (rev coded)","Well-Being index")
names(p.values_swb) <- names
adjusted.p_swb <- as.table(p.adjust(sort(p.values_swb),"BH"))
FDR_swb <- cbind(p.values_swb[names], adjusted.p_swb[names])
colnames(FDR_swb) <- c("Non-adjusted p-value","BH-adjusted p-value")
print(FDR_swb, digits = 4)  # columns 3 and 4
# For Table A10, treatment effect estimates and standard errors extracted from Table A7

# ------------- Table A11 [outgroup regard] ------------------
p.values_outgroup <- as.table(c(
  reg_res_mod3("ft",data)$p.value,
  reg_res_mod3("sd",data)$p.value,
  reg_res_mod3("cooperation",data)$p.value,
  reg_res_mod3("inchrct_data",data)$p.value,
  reg_res_mod3("othchrct_data",data)$p.value,
  reg_res_mod3("outgroup_index",data)$p.value,
  reg_res_mod3("outgroup_princomp",data)$p.value))
names<-c("feeling thermometer", "sd","cooperation","perception of out-group evaluations",
         "out-group traits","out-group regard (z-scores)","out-group regard (pc)")
names(p.values_outgroup) <- names
adjusted.p_outgroup <- as.table(p.adjust(sort(p.values_outgroup),"BH"))
FDR_outgroup <- cbind(p.values_outgroup[names], adjusted.p_outgroup[names])
colnames(FDR_outgroup) <- c("Non-adjusted p-value","BH-adjusted p-value")
print(FDR_outgroup, digits = 4)  # columns 3 and 4


# --------------------------------------
#     Robustness to outliers
# --------------------------------------

# We now examine the robustness of our findings to outliers by using Cook’s Distance:
# 1. first, we identify cases which are four times the mean value of Cook’s distance for all observations
# 2. then we exclude those cases, and re-estimate the model

var <- c("outgroup_index","outgroup_princomp","c_news","swb") # main outcome variables
sapply(var, function (i){
  sec <- lm(as.formula(paste0(i,"~","treat+as.factor(year)")), data=data)
  cookoutgroup_index<-cooks.distance(sec) 
  influential <- as.numeric(na.omit(names(cookoutgroup_index)[(cookoutgroup_index > 4*mean(cookoutgroup_index, na.rm=T))]))
  data$rownumber <- as.numeric(rownames(data))
  rep_rmalloutlier <- data[ ! data$rownumber %in% influential,]
  print(reg_res_mod3(paste0(i), rep_rmalloutlier))
})
# Results for third and fourth column (first two extracted from the results table)


#  ---------  MODERATION BY OFFLINE NETWORK

# ----------------------------------------------
#     FIGURE A2 (histogram of offline network)
# ----------------------------------------------

# Creating a simple sum score for offline network diversity, 
# for easier visualization and comprehension of the distribution of the offline network homogeneity
par(mar=c(5, 5, 5, 2))
data$offline_simple_sum <- NA
data[data$year=="21",]$offline_simple_sum <- data[data$year=="21",]$offline_network + data[data$year=="21",]$offline_exposure + data[data$year=="21",]$offline_active

c1 <- rgb(173,216,230,max = 255, alpha = 70, names = "lt.blue")
c2 <- rgb(255,192,203, max = 255, alpha = 80, names = "lt.pink")
a <- hist(data[data$ethn=="GR" & data$year=="21",]$offline_simple_sum, col="lightblue", xlim = c(0,16), freq = F, main="")
b <- hist(data[data$ethn=="TR" & data$year=="21",]$offline_simple_sum, col="lightpink", xlim = c(0,16), freq=F, main="")
plot(b,col=c1, xlim = c(0,16),ylim=c(0,1), freq=F,xlab="Offline network homogeneity (sum score)",main="")
plot(a, col=c2, xlim = c(0,16), ylim=c(0,1), add=TRUE,xaxt = 'n', yaxt = 'n', freq=F, box.lwd = 0, ,main = "Offline Network")
legend("topright", c("Turkish-speaking Cypriots", "Greek-speaking Cypriots"), col=c(c1, c2), lwd=8,cex=0.8, box.lwd = 0, bg = "white")


# ----------------------------------------------
#               TABLE A12
# ----------------------------------------------

# Column 1 of Table 12: Testing the interaction effect with the index of offline homogenous network (formed as a sum score)
interaction <- lm(outgroup_index ~ treat + offline_network_diversity + treat:offline_network_diversity, data[data$year=="21",])
summary(interaction)

# Column 2 of Table 12: Testing the interaction effect with the index of offline homogenous network (formed as a principal component score)
interaction_pc <- lm(outgroup_index ~ treat + offline_network_diversity_pc   + treat:offline_network_diversity_pc, data[data$year=="21",])
summary(interaction_pc)

# ---------------------------------------------------------------------------------------------------------
#     TABLE A13 (Moderation by offline network: Capital indicator as a proxy for offline heterogeneity)
# --------------------------------------------------------------------------------------------------------
data$capital <- ifelse(data$region=="nicosia" | data$region=="nicosia.tr", 0, 1)

# Interaction effect on the pooled data
capital <- lm(outgroup_index ~ treat + capital + treat:capital, data)
capital_coef <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,1]/sd(data[data$treat=="0" ,]$outgroup_index,na.rm=T)
capital_se <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_p <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,4]

capital_int <- lm(outgroup_index ~ treat + capital + treat:capital, data)
capital_int_coef <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,1]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_se <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_p <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,4]

capital_interactions <- data.frame(capital_coef, capital_se, capital_p, capital_int_coef, capital_int_se,capital_int_p)
capital_interactions


# Interaction effect on the 2020 data
capital <- lm(outgroup_index ~ treat + capital + treat:capital, data[data$year=="20",])
capital_coef <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,1]/sd(data[data$treat=="0" ,]$outgroup_index,na.rm=T)
capital_se <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_p <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,4]

capital_int <- lm(outgroup_index ~ treat + capital + treat:capital, data[data$year=="20",])
capital_int_coef <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,1]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_se <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_p <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,4]

capital_interactions <- data.frame(capital_coef, capital_se, capital_p, capital_int_coef, capital_int_se,capital_int_p)
capital_interactions

# Interaction effect on the 2021 data
capital <- lm(outgroup_index ~ treat + capital + treat:capital, data[data$year=="21",])
capital_coef <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,1]/sd(data[data$treat=="0" ,]$outgroup_index,na.rm=T)
capital_se <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_p <- coeftest(capital, vcov = vcovHC(capital, type = "HC1"))[2,4]

capital_int <- lm(outgroup_index ~ treat + capital + treat:capital, data[data$year=="21",])
capital_int_coef <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,1]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_se <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,2]/sd(data[data$treat=="0",]$outgroup_index,na.rm=T)
capital_int_p <- coeftest(capital_int, vcov = vcovHC(capital, type = "HC1"))[4,4]

capital_interactions <- data.frame(capital_coef, capital_se, capital_p, capital_int_coef, capital_int_se,capital_int_p)
capital_interactions


# ---------------------
#        FIGURE A3
# ---------------------
# Offline network diversity
data$offline_network_diversity <- NA
data[data$year=="21",]$offline_network_diversity <- scale(data[data$year=="21",]$offline_network) + scale(data[data$year=="21",]$offline_exposure) + scale(data[data$year=="21",]$offline_active)

quantile(data$offline_network_diversity, probs = c(0.25,0.50,0.75,1), na.rm=TRUE)

# Now we'll look at treatment effects within different quartiles of offline network homogeneity
reg_res_mod1y("outgroup_index",data=data[data$offline_network_diversity>1.86 & data$offline_network_diversity<=2.61 &data$year=="21",])
reg_res_mod1y("outgroup_index",data=data[data$offline_network_diversity>0.60 & data$offline_network_diversity<=1.86 &data$year=="21",])
reg_res_mod1y("outgroup_index",data=data[data$offline_network_diversity<=0.60 & data$offline_network_diversity>-1.40 &data$year=="21",])
reg_res_mod1y("outgroup_index",data=data[ data$offline_network_diversity<=-1.40 &data$year=="21",])

# Plotting
t <- tibble(x = 1:4, y = 1, z = 2, t=3)
colnames(t) <- c("estimate","sd","p.value","tau")
t$estimate <-c(-0.396, -0.205, 0.118 , 0.240)
t$sd<- c(0.313,0.312,  0.233 , 0.324)
t$tau <-c(0.25, 0.5, 0.75, 1)

ggplot(t, aes(x=tau, y=estimate)) + 
  geom_errorbar(aes(ymin=estimate-qnorm(.95)*sd, ymax=estimate+qnorm(.95)*sd, color=tau), width=.1) +
  geom_line() + geom_point()+
  theme(axis.text.x = element_text(size = 15), axis.title.x=element_text(size = 15),   axis.text.y = element_text(size = 15)) +
  xlab("Offline network heterogeneity") + ylab("Treatment effect")+ theme_light()+ theme(legend.position = "none")+ scale_x_continuous(breaks=c(0.25,0.5,0.75,1),labels = c("0%-25%","25%-50%","50%-75%","75%-100%"))+ geom_hline(yintercept=0, linetype="dashed", color = "red")

# Figure A4 was created using the nationally representative survey shared with the authors by the Centre for Sustainable Peace and Democratic Development.
# Access to that dataset cannot be granted without their permission, so please consider reaching out to either us or their team directly.

# ------------------------------------------------------
#      ADDITIONAL RESULTS [Appendix: Section G]
# ------------------------------------------------------
# --------- Victimization [Figure A4: left plot]

coef.vec1<- c(reg_res_mod2y("victim_narratives_1_downplay_ingroup",data[data$year=="21",])["estimate2"], reg_res_mod2y("victim_narratives_2_downplaying_outgroup",data[data$year=="21",])["estimate2"],
              reg_res_mod2y("victim_narratives_3_acknowledging",data[data$year=="21",])["estimate2"], reg_res_mod2y("victim_narratives_4_acknowledging",data[data$year=="21",])["estimate2"])
se.vec1<- c(reg_res_mod2y("victim_narratives_1_downplay_ingroup",data[data$year=="21",])["std.error2"], reg_res_mod2y("victim_narratives_2_downplaying_outgroup",data[data$year=="21",])["std.error2"],
            reg_res_mod2y("victim_narratives_3_acknowledging",data[data$year=="21",])["std.error2"], reg_res_mod2y("victim_narratives_4_acknowledging",data[data$year=="21",])["std.error2"])
par(mar=c(5, 15, 1, 1))
adjust = 0.10
layout(matrix(c(2,1),1,2.5), #in order to add variable categories and braces to left side of plot, 
       widths = c(2, 5.8))
var.names <- c("Downplaying \ningroup suffering","Downplaying \noutgroup suffering","Acknowledging \ningroup suffering","Acknowledging \noutgroup suffering")
y.axis <- c(length(var.names):1) 
unique(data$pol_orient)
#par(mar=c(5,1,5,1))
plot(coef.vec1, y.axis, type = "p", axes = F, xlab = "Standardized effect of deactivation", ylab = "", pch = 17,cex = 1,#plot coefficients as points, turning off axes and labels.
     xlim = c(-1,1), xaxs = "r", main = "Exposure to narratives..", col="black", cex.lab=1,mgp = c(3,1,0.2), cex.main=1.35)
segments(unlist(as.numeric(coef.vec1))-qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, unlist(as.numeric(coef.vec1))+qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, lwd = 1, col = c("black"))
axis(1, at = seq(-1,1,by=0.5), labels =c(-1,-0.5, 0,0.5,1), tick = T,#draw x-axis and labels with tick marks
     cex.axis = 1, mgp = c(4,2,1))#blackuce label size, moves labels closer to tick marks; zadnji u mgp spusta x liniju nize
axis(2, at = y.axis, label = var.names, las = 1, tick = T, mgp = c(1.6,.7,0),cex.axis = 1.15) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
abline(h = 0, v=0, col="gray",lty=2,lwd=2.5)

# [Figure A4: right plot]

coef.vec1<- c(reg_res_mod2y("victimization_1_othergroupssimilar",data)["estimate2"], reg_res_mod2y("victimization_2_uniquesuffering",data)["estimate2"],
              reg_res_mod2y("victimization_3_dontthinkmuch",data)["estimate2"])
se.vec1<- c(reg_res_mod2y("victimization_1_othergroupssimilar",data)["std.error2"], reg_res_mod2y("victimization_2_uniquesuffering",data)["std.error2"],
            reg_res_mod2y("victimization_3_dontthinkmuch", data)["std.error2"])
hist(data$sentiment_online)
par(mar=c(5, 15, 1, 1))
adjust = 0.10
layout(matrix(c(2,1),1,2.5), #in order to add variable categories and braces to left side of plot, 
       widths = c(2, 5.8))
var.names <- c("Inclusive \nvictimhood", "Competitive \nvictimhood","I do not think \nabout this")
y.axis <- c(length(var.names):1) 
plot(coef.vec1, y.axis, type = "p", axes = F, xlab = "Standardized effect of deactivation", ylab = "", pch = 17,cex = 1.35,#plot coefficients as points, turning off axes and labels.
     xlim = c(-1,1), xaxs = "r", main = "Group victimhood", col="black", cex.lab=1.15, cex.axis=1.25,mgp = c(3,1,0.2), cex.main=1.35)
segments(unlist(as.numeric(coef.vec1))-qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, unlist(as.numeric(coef.vec1))+qnorm(.975)*unlist(as.numeric(se.vec1)), y.axis, lwd = 1, col = c("black"))
axis(1, at = seq(-1,1,by=0.5), labels =c(-1,-0.5, 0,0.5,1), tick = T,#draw x-axis and labels with tick marks
     cex.axis = 1, mgp = c(4,2,1))#blackuce label size, moves labels closer to tick marks; zadnji u mgp spusta x liniju nize
axis(2, at = y.axis, label = var.names, las = 1, tick = T, mgp = c(1.6,.7,0),cex.axis = 1.35) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
abline(h = 0, v=0, col="gray",lty=2,lwd=2.5)